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Abstract 

Solving the Boltzmann - BGK equation, we investigate a flow generated by an 
infinite plate oscillating with frequency to. Geometrical simplicity of the prob- 
lem allows a solution in the entire range of dimensionless frequency variation 
< lot < oo, where r is a properly defined relaxation time. A transition from 
viscoelastic behavior of Newtonian fluid {lot — ► 0) to purely elastic dynamics 
in the limit lot — > oo is discovered. The relation of the derived solutions to 
microfluidics (high-frequency micro-resonators) is demonstrated on an example 
of a "plane oscillator". 



Introduction. During last two centuries, Newtonian fluid approximation was remark- 
ably successful in explaining a wide variety of natural phenomena ranging from flows 
in pipes, channels and boundary layers to the recently discovered processes in meteo- 
rology, aerodynamics , MHD and cosmology. With advent of powerful computers and 
development of effective numerical methods, Newtonian hydrodynamics remains at 
the foundation of various design tools widely used in mechanical and civil engineering. 
Since technology of the past mainly dealt with large ( macroscopic) systems varying 
on the length/ time-scales L and T, the Newtonian fluid approximation, typically 
denned by the smaliness of Knudsen and Weisenberg numbers Kn = X/L <C 1 and 
Wi = t/T « ±1 = KnMa < 1, was accurate enough. The length and time-scales A 
and r ~ A/c, are the mean -free path and relaxation time, respectively, . 
With recent rapid developments in nanotechnology and bioengineering, quantitative 
description of high-frequency oscillating microflows, i. e. flows where Newtonian 
approximation breaks down, became an important and urgent task from both basic 
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and applied science viewpoints. Modern micro (nano)- electromechanical devices 
(MEMS and NEMS), operating in the high- frequency range up to w = O(10 8 — 
10 9 Hz) (Wi = uot > 1), can lead to small biological mass detection, subatomic 
microscopy, viscometry and other applications [l]-[3]. The manufactured devices are 
so small and sensitive that adsorption of even tiny particles on their surfaces leads 
to detectable response in the resonator frequency peak (shift) and quality (Q) factor, 
thus enabling these revolutionary applications. Since both frequency shift and width 
of the resonance peak, depend upon properties of the resonator -generated flows, the 
microresonators may serve as sensors enabling accurate investigation of fundamental 
processes in microflows. 

Often, oscillating flows are a source of new and unexpected phenomena. For example, 
in a recent study of the high-frequency elctro-magnetic - field - driven nano -resonators 
(h x w x L ss 0.2 x x0.7 x 10/xm) Ekinci and Karabacak , observed a transition in the 
frequency dependence of inverse quality factor from 1/Q m oc 1/y/ui, expected 
in the hydrodynamic limit, to 1/Q oc 1/ou in the high-frequency (kinetic) limit [3]. 
Since parameter 7 oc W is proportional to the energy dissipation rate (W) into a 
surrounding gas, the discovered effect points to the frequency- independent dissipation 
rate W. No quantitative theory describing this transition has been developed. 
Recently, Shan et. al. performed a detailed numerical investigation of evolution of 
the initially prepared viscous shear layer defined by a unidimensional velocity field 
u(y, 0) = Mo cos(/co2/)i, where i is the unit vector in the x-direction. A novel transition 
from purely viscous decay to the viscoelastic dynamics was discovered for Tvk\ > 1 



In this paper we investigate deviations from Newtonian hydrodynamics, on an ex- 
ample of classic " Stokes' Second " problem of a flow generated by an infinite solid 
plate oscillating along rr-axis with velocity u(y — 0, t) — U cos cut [5]. First we use the 
hydrodynamic approximation, derived from kinetic equation by Chen et. al..[6] in the 
limit our <C 1. Then, derivation of equations valid in an infinite range < ujt < 00 
is presented and it is shown that in the high-frequency limit ujt — > 00, the dissi- 
pation rate W const. To apply the derived solution to experimental data on 
nano-resonators, we introduce a model system of plane oscillator and demonstrate 
the transition in the frequency dependence of the quality factor, in a quantitative 
agreement with experimental data of in Refs. [4] , [5] . 

The Boltzmann-BGK equation. Interested in the rapidly oscillating flows , where 
the Navier-Stokes equations break down, we consider the kinetic equation in the 
relaxation time approximation (RTA): 



In accord with Boltzmann's H-theorem , the initially non-equilibrium gas must mono- 



[4]- 
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tonically relax to thermodynamic equilibrium. This leads to the relaxation time an- 
zatz, qualitatively satisfying this requirement: 

J (2) 



t{v) 

In thermodynamic equilibrium the left side of kinetic equation is equal to zero and 

E(p,r) 

the remaining equation C = has a solution: f eq = constant x e e with the 
temperature 9 = const. In the mean -field approximation, valid close to equilibrium 
where all gradients are small, one replaces the variable v by v oc \JT jm. Given the 
mean-free path A, we have an estimate r = X/v ~ A/c = const where c is the speed 
of sound. In general, the relaxation time can be a non-trivial function of velocity 
gradients, position in the flow, external fields etc [7]. The equation (l)-(2) is the 
celebrated Boltzmann BGK equation widely used for both theoretical and numerical 
(Lattice Boltzmann Method) studies of non-equilibrium fluids [7]. This equation, 
with 

r = ^exp(-^£) (3) 

often considered as a generic equation of fluid dynamics, is valid when velocity gradi- 
ents are not large. It will be shown below that, due to a particular geometry studied 
in this paper, the largest magnitude of velocity gradient \d y u(0, t)\ < U/S with 5 — > oo 
when lot — > oo. Therefore, when the amplitude U is small enough, the typical gradi- 
ent is small. Moreover, this feature is responsible for the purely elastic response and 
total disappearence of viscosity from the problem. In turbulence theory this effect is 
called rapid distortioin (RD) limit. 

Chapman- Enskog expansion of BGK. Chen et.al (2003). Multiplying (1),(2) by v and 
integrating over v gives: 

dui dui Id . . 

^ + M ^ + p^' = (4) 
where the stress tensor, written for % ^ j is: 



a 



Ui)(vj - Uj) = / dv(vi - Ui)(vj - Uj)f(v,yi,t) (5) 



Usually, evaluation of the stress tensor is a difficult task. The simplified equation 
(l)-(2) for a single-particle distribution function allows calculation of nonlinear con- 
tributions to the momentum stress tensor = (vi — Ui)(vj — Uj) — \v 2 5ij. In a 
remarkable paper, Chen et. al. [6], based on thequation (l)-(2), formulated the 
Chapman-Enskog expansion in powers on dimensionless relaxation time e = rd Xi Uj. 
In this formulation: 
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V = eVi (6) 

9 d 2 d , d 
at at oti ot 2 
and the pdf is expanded in powers of e as: 

/ = / (0) + + e / (2) + • • • (8) 

Substituting these expressions into (1),(2) and equating the terms of the same order 
in e results in /(°) = f eq given by (3), and 

(A +v . v)/( o,_^ (9) 

( ^ + v .V)/« + ^ = -^ (10) 
etc. The mean of any fluctuating variable is then calculated as : 

~My) = A (0) + A w + A (2) + • • • (11) 



where 



AW = J /< n U(v)v 



To illustrate the main results and simplify notation in what follows we set the tem- 
perature 9 = const. The calculation of Chen (2003) gives: 



« ~f i0 SiA(vi - - Uj) - (12) 



and 



f (2) = _2 r 2 / (0) [K - «<)^:(5ii - • u<y + 0(( Up - Mp )K - u,)^*)] (13) 

It is important to stress that the last contribution to the right side of (13) involves 
even powers of v — u and various products of SV,-. The result is (Chen (2003)): 



an = 2T9Sij -2r9(dt + u- V)(r£, 



— 4r 2 9[S ik Skj — ^SijSkiSki] + 2T 2 9[S ik VL kj + S jk Q k i] (14) 



where the rate of strain = ^[f^ 1 + and the vorticity tensor is defined as, 

o — ird"» _ ®m] 

The first term in the right side of (14), resulting from the first order Chapman-Enskog 
(CE) expansion, corresponds to the familiar Navier-Stokes equations for Newtonian 
fluid and the non-linear (non - Newtonian) corrections, given by the remaining terms, 
are generated in the next , second, order. The constitutive relation (14) is quite 
complex and, in general, can be attacked by numerical methods only. However, it is 
greatly simplified in an important class of simple unidirectional flows. 
Stokes' Second Problem is formulated as follows: The flow of a fluid filling half -space 
< y is generated by a solid plate at y = moving along the x-axis with velocity 
U p (0, t) = U cos out. Since velocity components in y and z-directions are equal to zero, 
we have to solve the equations for the x-component of velocity field u(y, t) only. Due 
to geometry of the problem 

u x = u — u(y, t); u y = v = 0; u • V = 0; d x = (15) 

Thus, since i ^ j, we are interested in a x>y . In this case, the equation of motion 
corresponding to the stress tensor (14) is very simple: 



Kl-^fe (16) 



du d d 2 u 

dt = ^ ~ T dt'df 

In the limit out — > 0, this equation is to be solved subject to the no-slip boundary 
condition, u(0,t) = U cos out; lim^oo u(y, t) = 0. Seeking a solution, satisfying the 
boundary condition at y — > oo, as: u = TZ{(j){ii)e~' lwt ) gives: = Be~~s and in the low 
frequency limit Wi — tuj <^ 1: 

u(y, t) = Uexp{-^{1 - T -^)y) cos{out - ^(1 + ^-)y) (17) 



In a classic case out = 0, the flow is characterized by a single scale ^2v/ou describ- 
ing both "penetration depth" and wave-length of transverse waves radiated by the 
oscillating plate. We can see that the non-newtonian O(rou) contribution leads to 
formation of two different length-scales: the 0{rou) increasing penetration depth and 
by the wave-length decreasing by the same magnitude. This is a qualitatively new 
feature of a non-newtonian flow. Now we calculate the dissipation rate. The force 
acting on the unit area of the wall at y = : 

8 .du(0,t) TT , r , 7T S TUJ , 7T N , 

^ ~ T di ) ^T = "^V^WKM + i ) + -y cos M - 4 )] ( 18 ) 

is phase-shifted relative to velocity field u(0, t) = U cos out. This result differs from 
its classic Newtonian counterpart by an 0(tou/2) shift. The mean energy dissipated 
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per unit time per unit are of the plate is calculated if we multiply (18) by —U cos cut 
and integrate over one cycle with the result: 



W^) = ^f(l + ^W(0,aO (19) 

In a simple gas close to thermodynamic equilibrium, the mean-free path A ~ ^ L Sn 

where £ oc d 2 is a scattering cross-section of the molecules and d is a typical 
scale of intermolecular interaction. 

Large deviations from Newtonian fluid mechanics. Here we somewhat modify the 
procedure developed by Chen (2003). It is clear that as r — > 0, the expansion gives the 
classic Stokes results for the flow of oscillating plate. Moreover, in the limit too <^1, 
the second -order in S correction to (14) disappears due to the symmetries of the 
problem defined by (7). The expression (14) includes a time derivative, which for the 
oscillating flow we are interested in this paper, introduces an additional dimensionless 
parameter Wi — tuj into expansion. Therefore, the CE expansion is in fact expansion 
in two dimensionless parameters Kn = rd y u and Wi. As will be shown below, in the 
simple-geometry oscillating flows, these parameters are quite different : as tuj — > oo, 
the second parameter Kn — > 0. Thus, while neglecting the small, 0(V 2ri w) with n > 1, 
"Burnett contributions", we will attempt to sum up the entire series in Wi = tuj . 
The Boltzmann equation is: 

^l + Va ^L + L = ll (20) 

dt a dx a T T 

Multiplying (20) by (i>j — Ui)(vj — Uj) and integrating over v, we derive the equation 
for the stress tensor ( % ^ j): : 

i)cJ ''> + a JL = _I f dv ( v . _ u .){y. _ Uj )v a ^- (21) 



dt t p J ; dx a 

or 



do~ij <7ij duj d — — 

+ u • Va, tj + — = -2a iya — - - «a)(u< - Ui)\Vj - Uj) - o-ijV ■ u (22) 

Due to the basic symmetries of the problem (7), this equation is simplified: 

13 + — = -2<7j iQ |^- - -£—(v a - u a )(vi - Ui)(vj - uj) (23) 



dt t ' dx n dx, 



a 



In the unidirectional flow, we are interested in only d y u ^ and therefore, 
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°i, a d Xa U = CTyyOyU (24) 

The remaining term in the right side of equation (15) can also be simplified leading 



to: 



dcr y , x . °y,x ~~du d 



at + T^- 2 %-^-*> < 25 > 



where = o y<lj — Vy — — ^ . This equation is formally exact. 
Our goal now is to evaluate 



Eij = ^~ (Uq " Ua){Vt ~ Ui)iVj ~ Uj) = E ^° } + ^ + ^ + " ' (26) 



This can be done readily with relations (3), (11)- (13) derived in Chen (2003). We see 
if = and Eg> 



that iff = and Y$ = 0. Evaluated on from (13) 



(27) 

where A\ = J f°(v)A(v)dv. It will be shown a'posteriory that in both limits tuj — > 
and rc<j — > oo, the second order contribution 



r^-< 



6d T u 5 2 



(28) 



where 5 is the width of the viscous layer. 

As we saw, the high powers of S disappear from the second-order relation (14). Let 
us see that this is a general phenomenon valid to all orders. Let us consider a general 
n th -order contribution to the stress tensor a x>y \ 

du ai du a2 du a3 du x <9w Qn _ 1 du an 
dx a:i dxy dx anl 9x an dx a2 dx ai 

where the summation is carried out over the randomly distributed Greek subscripts. 
By the symmetry 

du x du a du dv . . 

= (30) 



dx an dx ai dy dx ai 

because v — 0. Denoting a XjV = a we, based on the above considerations obtain 
equation valid in both low and large-frequency limits: 

du da da du 

m=-dy- ; T m +a = - u dy- (31) 
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Comparing this result with the equation (16), we conclude that constitutive equa- 
tion (31) is a resummation of the Chap man- Enskog expansion used in the low-order 
derivation of (14). Indeed, the Fourier-transform of the second equation in (31) is 
equal to (1 + icuT)a(y,cu) = which in the limit cut — > 0, coinsides with (16). The 
equations (31) give: 

d 2 u du d 2 u . . 

T w + m = V W (32) 

Boundary conditions. Interested in the limit cut > 1, we have to be careful with 
the choice of boundary conditions (BC). The BC used in this work is u(0,t) = 
U{ut) coscut = ol(u)t)U coscut. The "slip factor < a(wr) < 1 has been recently 
investigated in detailed numerical simulations where is was shown that for cut < 1, 
a(Lur) ~ 1 and it rapidly decreases to o.(cut) ks 0.3 — 0.2 for 1 < cut < 100 [9]. Deal- 
ing with the linear equation, to remove the slip factor from consideration, we can, in 
principle, introduce the normalized velocity field u/U(cut), solve equations and then 
recover u(y,t). Below, we use the no-slip boundary condition. 

solution. The solution to equation (32) is found readily: u = 7Z[e~ tuJt e~~s] with 



i 



where 



)Jj-Vl — an d: 



u = Ue cos(cjt - j-) (33) 



1/4 = (l+cu 2 r 2 )3^[ C os(itan- 1 (^r))±sin(itan- 1 (cur))] (34) 

In the limit wr — > 0, this solution tends to the expression (17) derived above. As we 
see, the non-Newtonian Second Stokes problem can be characterized by two different 
length-scales, penetration length 5_ and the wave-length 5 + . In the limit cut — > oo, 
the penetration length tends to infinity and the dominating dissipation mechanism 
is the wave generation. (This limits the velocity gradients in this problem). To 
calculate the force acting on the unit area of the plate, we notice that in accord with 
the differential equations (31) 

<7(0,*)=pi/exp(--) /* ^ ( ° ,A) exp(-)d- (35) 
r J~oo ay r t 

where z/^i = u(-^ + siff ) and 

cr(0, t) = pv- — ^ [— — (cosa;t + cut sin cut) + — (—cut cos cut + smut)] (36) 

1 + CU 2 T 2 0_ 6+ 
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Figure 1: Dissipation rate 2W(t,lu) as a function of u (expression (16); (arbitrary 
units). 

where p is the density of a fluid. To obtain the dissipation rate per unit time per unit 
area of the plate, we calculate W = —u(0,t)a(0,t) averaged over a single cycle. The 
result is: 



W(t,u) 



1 pU 2 



v cL + 6+' 



(37) 



2 1 + UJ 2 T 2 

The plot of the normalized dissipation rate as a function of ujt is shown on Fig.l . 
As ujt — * 0, we recover the previous result W/y/u oc (1 + ^). In the opposite limit 
ujt — > oo, the saturation of the curve W{ujt) is predicted. In this range, the kinetic 
energy of the plate oscillations is mainly dissipated into the wave radiation. 



"Plane oscillator" . Consider a massless spring of stiffness k with two infinite plates 
of height h attached to it. The losses in a spring are neglected and the friction force, 
acting on the plates, is the only source of energy dissipation. The equation of motion 
of this "plane oscillator" driven by a force R(t) is: 

x a + ix t + kx = R(t) (38) 

where jx t = 7ii(0, t) = 2pa(0,t)/m s is the friction force acting on a unit mass 
(m s = p p h ) of a plate surface and the factor 2 accounts for the force acting on both 
(top and bottom) surfaces of the plate. If the resonance I(uo) is sharp enough, so 
that 5uj/ujq <^ 1, then we can calculate the damping (friction) force acting on the 
oscillator using the theory developed above. In Newtonian limit (hereafter, we omit 

the subscript ) ujt — > 0, the expression (36) gives 7 = ^ p:puj / {2m 2 s ) oc ^Joj and the 
quality factor (the definition is given below) Q ~ w/7 oc y/cu. In the opposite (kinetic) 
limit uj — > 00 andr = const: 7 — > pujT J ^ / (m s {l+uj 2 T 2 ) 4) oc uj° = const = 0(pc/m s ) 
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This relation can be understood as follows. We are interested in a somewhat un- 
usual case with the rapidly varying deterministic velocity field u(0, t) and the slow - 
varying ("frozen") molecular motion. Thus, the long -time-averaging, giving u — 0, 
is unphysical and the stress on a plane y = is calculated by averaging over the 
plate surface: a xy = (u(0,t) + c) x (u(0,t) + c) y ~ u(0, t)c s , where c s is the mean 
velocity (y-component) of molecules colliding with solid surface. In this limit, the 
tangential component of the force acting on a unit surface of a solid (per/unit mass) 
is: F t = ^u(0,t) and 7 = The former relation, derived here for a rapidly 

oscillating plate, is similar to the one appearing in a problem of a piston instan- 
taneously accelerated in an ideal gas. Thus, in this limit the inverse quality fac- 
tor 1/Q = , 7 ps 7/aj oc 1/uj. For a gas of a given temperature p oc p and 

c oc \JO/m mo i = const, we have ^ oc p/uj 

To obtain the relation valid in the entire range of frequency variation we have from 
(31): 



T = 2 ~T77 ^TaK 1 + ut p) cost-tan' 1 u;t p ) - (1 - ut p ) sin(-icm~W p )] (39) 

Pph(l + UJ Z T^)i * * 

Our goal now is to find the relaxation time r p responsible for relaxation to equilib- 
rium in the immediate vicinity of the rapidly oscillating solid plate. In a standard 
equilibrium situation, when A and r ps A/c 1/uj are the smallest length and 
time-scales in the system, the relaxation times in the bulk and in the immediate 
vicinity of the surface are more or less equal. In a rapidly oscillating flow, we are 
interested in, this is not so. Indeed, even in the air at normal conditions, the bulk 
t ~ A/c Ps 10~ 9 sec and in the modern microresonators, where uo ~ 10 9 sec -1 , the 
inequality ujt <^ 1 is hardly satisfied. Moreover, in the low pressure devices where 
1CT 3 < p < 10 3 torr , the bulk relaxation time is huge, so that r p must be calculated 
from a theory taking into account strongly non-equilibrium effects. In the absence of 
such theory, we use the results of experimental data on a driven microbeam obtained 
by Ekinci and Karabacak [3] who covered an extremely wide range of parameter 
variation: 10° < p < lOHorr, 10 6 < u < 0.6 • 10 9 sec -1 , h = 2 • 10~ 5 cm. Under 
the normal conditions for u = 0.6 • 10 9 sec _1 , the observed damping parameter was 
7 ps 1.5 • 10 6 sec _1 . (This relation is the result of the measured 1/Q = 7 /uo ps 3 ■ 10~ 3 
and c^o = 0.6 • 10 9 sec~ 1 .) Denoting x = ujt p and introducing a geometric factor g 
accounting for the geometry differences between theoretical and experimental setups, 
we have: 

X 1 

u l + x) cos(-tan _1 x) — (1 — x) sin(-tan~ 1 x)] = 1 (40) 



4(l + x 2 )i v 2 v 2 
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In this paper we considered a case of the simplest possible geometry. In the nano- 
technological applications the resonator amplitude is 0(lCT 6 cm) is much smallest 
than the smallest linear dimension of the beam (cantelever). Thus, the results pre- 
sented here may be not too geometry -sensitive and can be readily generalized to 
different geometries. For example, the solution to the problem of an infinite cylinder 
oscillating along its axis can be expressed in terms of the Bessel functions leading to 
0(1) variations of the coefficient in front of (40). Here, to qualitatively illustrate the 
origins of the geometric factor g, we use very simple considerations. As was stated 
above, mass per unit area of a plate of height h and length / is : m s = p p l 2 h/l 2 = p p h. 
A simple calculation of m s for circular cylinder gives: m s = pa/2, where a is the 
radius. In this case, g cy i ~ Ad p i ate . For a beam of height h and width w, we have 
m s = Pp 2(w+h) ■ Thus, the dissipation per unit area may be a relatively universal 
property [10]. The equation (19), written for a particular magnitude of frequency 
(u> = 0.6 • 10 9 sec _1 ) and other parameters of experiment of Refs.[3], [10]) has solu- 
tions for all g > 3.6. For example, g — 4 and g = 6.28 give x ~ 1.8 and x ~ 4.8 
corresponding to r pa 2nsec and r pa 8nsec, respectively. Choosing g — 4 we, using 
the experimental data for normal pressure p = 800torr as a calibration point take 
r p (800torr) pa 2nsec [3]. Then, based on kinetic theory, we substitute r p = ^^-nscc 
into (39) and obtain 7 and quality factor Q(u>,p) in a good agreement with the results 
of Ekinci and Karabacak [3] in a wide range of both frequency and pressure variation 
(10 6 < u < 0.6 x l^Hz; 100 < p < lOOOtorr. ) 

In a kinetic limit 7 = const, and introducing G(x) = [(1 + x) cos(| tan -1 x) — (1 — 
x) sin(i tan~ 1 x)]/ (1+x 2 ) 4, the relation (39) give: a^fu:G{x) = 1, where the coefficient 
a = a(p, 9, m s , p, g) = const includes all parameters of the system. It is easy to check, 
that solution to this equation x oc cu defines an asymptotic kinetic time t p = const 7^ r 
characterizing relaxation to equilibrium in the high-frequency limit wr — > 00. Thus, 
the equilibrium and non-equilibrium (high-frequency) relaxation times may differ by 
a constant 0(1) factor only. The inverse quality factor \/Q{uj) ps 7/u;, given by (39), 
is shown on Fig. 2. The most detailed comparison of the results based on relation 
(39) with experimental data is presented in Ref. 10. 

In a dense liquid where ojt 1, the relation (39) gives: 



l J P2(w+h) 



and in water [y ps 0.01^-, p/p p ^ 2, u pa 0.6 ■ 10 9 sec 1 ; h = 2 ■ 10 5 cm and 
w — 7 • 10~ 5 cm), corresponding to the experimental set up of Ref. 11, we predict 

I 2 2 

7 pa 1.6 • lOSecT 1 and Q pa ^— — — pa 4, in a good agreement with the data [11]. 
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10 5 10 6 10 7 10 8 10 ! 

Frequency a> (1/s) 



Figure 2: Inverse quality factor l/Q m 7/0* = j^j=G(x) ~ 60^ g§^G(x) (relation 
(18)) vs frequency u in the presure interval 100 < p < lOOOtor. Relaxation time 
r p (800torr) = 2nsec. For remaining parameters, see text. 

Summary and conclusions. In this paper a theory of a flow generated by the oscillating 
plate, valid in a wide range of variation of dimensionless frequency (Weissenberg 
number) < lot < 00, is presented. The solution to kinetic BGK equation, describing 
transition between viscoelastic and purely elastic regimes, has been found. The results 
obtained for a simple "plane oscillator", first introduced in this paper, were favorably 
compared with the outcome of experiments on nanoresonators of Refs. [3], [10] where 
this transtion has experimentally been observed. 

We would like to thank K. Ekinci and D. Karabacak for sharing with us their yet 
unpublished data on microresonators. We benefited from the most interesting and 
stimulating discussions with them and R. Benzi, H. Chen, R. Zhang, X. Shan, F. 
Alexander, S. Succi, I. Karlin and V. Steinberg. 
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